source("read.data.R")
source("util.R")
source("analy.model.R")
source("analy.model.selection.R")
source("analy.QRE.general.R")
source("analy.plot.R")

library(lme4)
library(lmerTest)

# e0<-par.setting[1] 
# ea<-par.setting[2]
# inv.cost<-par.setting[3]
# p0<-par.setting[4]
# p1<-par.setting[5]
# penalty<-par.setting[6]
# 
# gamma.pay<-par.behavior[1]
# gamma.r<-par.behavior[2]
# gamma.a<-par.behavior[3]
# gamma.i<-par.behavior[4]

time.start<-proc.time()

library(parallel)
cluster<-makeCluster(6)
clusterExport(cl=cluster,ls(),envir=.GlobalEnv)

cn2<-c("gamma pay","gamma ransom","gamma attack","gamma invest","follow me inv","pay social norm","pay fairness","inv social norm","follow you inv")
cn.list<-list(cn2)

treat<-"Base"
data<-read.file.index(1)
data<-subset(data,data[,"Period"]>1)
data.base<-subset(data,data[,"Treatment"]==treat) 

par.setting<-c(data.base[1,"Data1"],data.base[1,"Outside"],data.base[1,"Investment1"],0.8,0.5,data.base[1,"Penalty"],0,0,0)
par.b<-c()

par.b<-convert.par.stl(ms[[2]]$par,model.spec=ms[[1]])
par.b<-as.numeric(par.b)
par.b<-ifelse(is.na(par.b),ms[[7]],par.b)

par.b[5]<-0
par.b[9]<-0

# 
# r<-pmax(1,data.base[,"ransom"])
# inv<-ifelse(data.base[,"attack.decision"]==1,data.base[,"i_1"],data.base[,"i_2"])
# pay<-ifelse(data.base[,"attack.decision"]==1,data.base[,"p_1"],data.base[,"p_2"])
# 
# pp1<-prob.pay(c(1:100),1,par.setting=par.setting,par.behavior=par.b)
# pp0<-prob.pay(c(1:100),0,par.setting=par.setting,par.behavior=par.b)
# #pp.model<-mean(na.omit(pp[r]))
# #pp.fit<-mean(na.omit(c(data.base[,"p_1"],data.base[,"p_2"])))
# 
# data.pay<-cbind(r,inv,pay)
# data.pay<-subset(data.pay,(!is.na(data.pay[,"r"]))&(!is.na(data.pay[,"pay"])))
# 
# pp.data<-c(mean(data.pay[data.pay[,"inv"]==0,"pay"]),mean(data.pay[data.pay[,"inv"]==1,"pay"]))
# pp.model<-c(mean(pp0[data.pay[,"r"]][data.pay[,"inv"]==0]),mean(pp1[data.pay[,"r"]][data.pay[,"inv"]==1]))
# 
# pp.summary<-cbind(c(0,1),pp.model,pp.data)
# dimnames(pp.summary)[[2]]<-c("inv","model","data")
# 

qre.r.calc<-function(par,par.setting=NULL) {
  
  par.b<-par
  
  u.m1<-matrix(nrow=2,ncol=2)
  u.m2<-matrix(nrow=2,ncol=2)
  
  for(i in 0:1) for(j in 0:1) {
    u.m1[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
    u.m2[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
  }
  
  dimnames(u.m1)<-list(c("0","1"),c("0","1"))
  dimnames(u.m2)<-list(c("0","1"),c("0","1"))
  
  qre.i<-find.qre(c(par.b[4],par.b[4]),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
  # qre.i<-find.qre(c(0.01,0.01),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
  
  qi1<-qre.i[[1]]
  qi2<-qre.i[[2]]
  inv.summary<-c("no-no","no-yes","yes-no","yes-yes")
  inv.summary<-cbind(inv.summary,c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2]))
  inv.p<-c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2])
  
  
  qrr1<-qre.r(1,par.setting=par.setting,par.behavior=par.b)
  qrr0<-qre.r(0,par.setting=par.setting,par.behavior=par.b)
  
  qra00<-qre.att(c(0,0),par.setting=par.setting,par.behavior=par.b)
  qra01<-qre.att(c(0,1),par.setting=par.setting,par.behavior=par.b)
  qra10<-qre.att(c(1,0),par.setting=par.setting,par.behavior=par.b)
  qra11<-qre.att(c(1,1),par.setting=par.setting,par.behavior=par.b)
  
  prob.r0<-inv.p[1]*sum(qra00[2:3,2])+inv.p[2]*qra01[2,2]+inv.p[3]*qra10[3,2]
  prob.r1<-inv.p[2]*qra01[3,2]+inv.p[3]*qra10[2,2]+inv.p[4]*sum(qra11[2:3,2])
  
  prob.sum<-prob.r0+prob.r1
  
  prob.r0<-prob.r0/prob.sum
  prob.r1<-prob.r1/prob.sum
  
  er<-sum(prob.r0*qrr0[,1]*qrr0[,2]+prob.r1*qrr1[,1]*qrr1[,2])
  er   
} 



par.b.list<-list(par.b)
data.qre.r.inv<-calculate.multiple(qre.r.calc,par.b.list,key.index=8,range=c(0,50),n=20,par.setting=par.setting)
data.qre.r.pay<-calculate.multiple(qre.r.calc,par.b.list,key.index=6,range=c(0,50),n=20,par.setting=par.setting)

jpeg(file="output/sim.r.jpg")
plot.multiple(data.qre.r.inv[,1],cbind(data.qre.r.inv[,2],data.qre.r.pay[,2]),c("investment","payment refusal"),xlab="treatment utility",ylab="expected ransom")
dev.off()

#plot.multiple(data.qre.r[,1],data.qre.r[,c(2,3)],c("target did not invest","target invested"),xlab="not-pay social norm",ylab="expected ransom")

#par.b.list<-list(c(par.b,0),c(par.b,1))
#data.qre.r<-calculate.multiple(qre.r.calc,par.b.list,key.index=8,range=c(0,50),n=20,par.setting=par.setting)

#plot.multiple(data.qre.r[,1],data.qre.r[,c(2,3)],c("target did not invest","target invested"),xlab="not-pay social norm",ylab="expected ransom")


# 
# qrr<-cbind(qrr1[,1],rep(0,nrow(qrr1)),rep(0,nrow(qrr1)))
# 
# for(i in 1:nrow(data.base)) {
#   if(data.base[i,"attack.decision"]>0) {
#     r<-max(1,data.base[i,"ransom"])
#     
#     old.payment.outcome<-ifelse(is.na(data.base[i,"old_payment_outcome"]),0,data.base[i,"old_payment_outcome"])
#     par.setting[7]<-old.payment.outcome
#     
#     ii<-ifelse(data.base[i,"attack.decision"]==1,data.base[i,"i_1"],data.base[i,"i_2"])
#     qrr.i<-qre.r(ii,par.setting=par.setting,par.behavior=par.b)
#     
#     qrr[r,3]<-qrr[r,3]+1
#     if(r<1) browser()
#     if(r>100) browser()
#     
#     qrr[,2]<-qrr[,2]+qrr.i[,2]
#   }
# }
# 
# qrr[,2]<-qrr[,2]/sum(qrr[,2])
# qrr[,3]<-qrr[,3]/sum(qrr[,3])
# 
# dimnames(qrr)[[2]]<-c("r","model","data")
# 
# qrr.plot<-plot.smooth(qrr,n=10)
# 
# ff<-paste("output/figure.",format(Sys.time(), "%Y_%m_%d_%H_%M_%S"),".",treat,".",model.list[model.index],".jpg",sep="")
# 
# m.lab<-paste(treat," ",model.list[model.index],sep="")
# plot.multiple(qrr.plot[,1],qrr.plot[,c(2,3)],c("model","data"),file=ff,main=m.lab,xlab="ransom",ylab="QRE prob")
# 

# att.summary<-cbind(c(0,1,2),c(0,0,0),c(0,0,0))
# for(i in 1:nrow(data.base)) {
#    inv.list<-c(data.base[i,"i_1"],data.base[i,"i_2"])
#    qra<-qre.att(inv.list,par.setting=par.setting,par.behavior=par.b)
#    att.summary[,2]<-att.summary[,2]+qra[,2]
#    att.summary[(att.summary[,1]==data.base[i,"attack.decision"]),3]<-att.summary[(att.summary[,1]==data.base[i,"attack.decision"]),3]+1
# }
#  
# att.summary[,2]<-att.summary[,2]/sum(att.summary[,2])
# att.summary[,3]<-att.summary[,3]/sum(att.summary[,3])
#  
# dimnames(att.summary)[[2]]<-c("att","model","data")
# 
# 

qre.i.calc<-function(par,par.setting=NULL) {

  par.b<-par
  
  u.m1<-matrix(nrow=2,ncol=2)
  u.m2<-matrix(nrow=2,ncol=2)
  
  for(i in 0:1) for(j in 0:1) {
    u.m1[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
    u.m2[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
  }
  
  dimnames(u.m1)<-list(c("0","1"),c("0","1"))
  dimnames(u.m2)<-list(c("0","1"),c("0","1"))
  
  qre.i<-find.qre(c(par.b[4],par.b[4]),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
  # qre.i<-find.qre(c(0.01,0.01),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
  
  qi1<-qre.i[[1]]
  qi2<-qre.i[[2]]
  inv.summary<-c("no-no","no-yes","yes-no","yes-yes")
  inv.summary<-cbind(inv.summary,c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2]))
  
#  browser()
  
  (qi1[2,2]+qi2[2,2])/2
}

par.b.list<-list(c(par.b))
data.qre.i.inv<-calculate.multiple(qre.i.calc,par.b.list,key.index=8,range=c(0,50),n=20,par.setting=par.setting)
data.qre.i.pay<-calculate.multiple(qre.i.calc,par.b.list,key.index=6,range=c(0,50),n=20,par.setting=par.setting)

jpeg(file="output/sim.i.jpg")
plot.multiple(data.qre.i.inv[,1],cbind(data.qre.i.inv[,2],data.qre.i.pay[,2]),c("investment","payment refusal"),xlab="treatment utility",ylab="investment rate")
dev.off()

qre.pay.calc<-function(par,par.setting=NULL) {
  
  par.b<-par
  
  u.m1<-matrix(nrow=2,ncol=2)
  u.m2<-matrix(nrow=2,ncol=2)
  
  for(i in 0:1) for(j in 0:1) {
    u.m1[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
    u.m2[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
  }
  
  dimnames(u.m1)<-list(c("0","1"),c("0","1"))
  dimnames(u.m2)<-list(c("0","1"),c("0","1"))
  
  qre.i<-find.qre(c(par.b[4],par.b[4]),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
  # qre.i<-find.qre(c(0.01,0.01),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
  
  qi1<-qre.i[[1]]
  qi2<-qre.i[[2]]
  inv.summary<-c("no-no","no-yes","yes-no","yes-yes")
  inv.summary<-cbind(inv.summary,c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2]))
  inv.p<-c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2])
  
  
  qrr1<-qre.r(1,par.setting=par.setting,par.behavior=par.b)
  qrr0<-qre.r(0,par.setting=par.setting,par.behavior=par.b)
  
  qra00<-qre.att(c(0,0),par.setting=par.setting,par.behavior=par.b)
  qra01<-qre.att(c(0,1),par.setting=par.setting,par.behavior=par.b)
  qra10<-qre.att(c(1,0),par.setting=par.setting,par.behavior=par.b)
  qra11<-qre.att(c(1,1),par.setting=par.setting,par.behavior=par.b)
  
  prob.r0<-inv.p[1]*sum(qra00[2:3,2])+inv.p[2]*qra01[2,2]+inv.p[3]*qra10[3,2]
  prob.r1<-inv.p[2]*qra01[3,2]+inv.p[3]*qra10[2,2]+inv.p[4]*sum(qra11[2:3,2])
  
  prob.sum<-prob.r0+prob.r1
  
  prob.r0<-prob.r0/prob.sum
  prob.r1<-prob.r1/prob.sum
  
  pp1<-prob.pay(c(1:100),1,par.setting=par.setting,par.behavior=par.b)
  pp0<-prob.pay(c(1:100),0,par.setting=par.setting,par.behavior=par.b)
  
  pp<-sum(prob.r0*pp0*qrr0[,2]+prob.r1*pp1*qrr1[,2])
#  browser()
  pp  
} 

par.b.list<-list(c(par.b))
data.qre.pay.inv<-calculate.multiple(qre.pay.calc,par.b.list,key.index=8,range=c(0,50),n=20,par.setting=par.setting)
data.qre.pay.pay<-calculate.multiple(qre.pay.calc,par.b.list,key.index=6,range=c(0,50),n=20,par.setting=par.setting)

jpeg(file="output/sim.pay.jpg")
plot.multiple(data.qre.pay.inv[,1],cbind(data.qre.pay.inv[,2],data.qre.pay.pay[,2]),c("investment","payment refusal"),xlab="treatment utility",ylab="payment rate")
dev.off()


qre.attack.calc<-function(par,par.setting=NULL) {
  
  par.b<-par
  
  u.m1<-matrix(nrow=2,ncol=2)
  u.m2<-matrix(nrow=2,ncol=2)
  
  for(i in 0:1) for(j in 0:1) {
    u.m1[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
    u.m2[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
  }
  
  dimnames(u.m1)<-list(c("0","1"),c("0","1"))
  dimnames(u.m2)<-list(c("0","1"),c("0","1"))
  
  qre.i<-find.qre(c(par.b[4],par.b[4]),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
  # qre.i<-find.qre(c(0.01,0.01),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
  
  qi1<-qre.i[[1]]
  qi2<-qre.i[[2]]
  inv.summary<-c("no-no","no-yes","yes-no","yes-yes")
  inv.summary<-cbind(inv.summary,c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2]))
  inv.p<-c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2])
  
  
  qrr1<-qre.r(1,par.setting=par.setting,par.behavior=par.b)
  qrr0<-qre.r(0,par.setting=par.setting,par.behavior=par.b)
  
  qra00<-qre.att(c(0,0),par.setting=par.setting,par.behavior=par.b)
  qra01<-qre.att(c(0,1),par.setting=par.setting,par.behavior=par.b)
  qra10<-qre.att(c(1,0),par.setting=par.setting,par.behavior=par.b)
  qra11<-qre.att(c(1,1),par.setting=par.setting,par.behavior=par.b)
  
  prob.r0<-inv.p[1]*sum(qra00[2:3,2])+inv.p[2]*qra01[2,2]+inv.p[3]*qra10[3,2]
  prob.r1<-inv.p[2]*qra01[3,2]+inv.p[3]*qra10[2,2]+inv.p[4]*sum(qra11[2:3,2])
  
  prob.sum<-prob.r0+prob.r1
  
  prob.r0<-prob.r0/prob.sum
  prob.r1<-prob.r1/prob.sum
  
  pp1<-prob.pay(c(1:100),1,par.setting=par.setting,par.behavior=par.b)
  pp0<-prob.pay(c(1:100),0,par.setting=par.setting,par.behavior=par.b)
  
  # pp<-sum(prob.r0*pp0*qrr0[,2]+prob.r1*pp1*qrr1[,2])
  # browser()
  
  pp<-inv.p[1]*sum(qra00[2:3,2])+inv.p[2]*sum(qra01[2:3,2])+inv.p[3]*sum(qra10[2:3,2])+inv.p[4]*sum(qra11[2:3,2])
  pp  
} 

par.b.list<-list(c(par.b))
data.qre.attack.inv<-calculate.multiple(qre.attack.calc,par.b.list,key.index=8,range=c(0,50),n=20,par.setting=par.setting)
data.qre.attack.pay<-calculate.multiple(qre.attack.calc,par.b.list,key.index=6,range=c(0,50),n=20,par.setting=par.setting)

jpeg(file="output/sim.attack.jpg")
plot.multiple(data.qre.attack.inv[,1],cbind(data.qre.attack.inv[,2],data.qre.attack.pay[,2]),c("investment","payment refusal"),xlab="treatment utility",ylab="attack rate")
dev.off()


stopCluster(cluster)

time.end<-proc.time()
tt<-time.end[3]-time.start[3]
tt<-tt/60
cat("time elapsed = ",tt," minutes\n")
save.image()